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We investigate close binary neutron stars in quasiequilibrium states in a general relativistic frame- 
work. The configurations are numerically computed assuming (1) existence of a helicoidal Killing 
vector, (2) conformal flatness for spatial components of the metric, (3) irrotational velocity fleld 
for the neutron stars and (4) masses of neutron stars to be identical. We adopt the polytropic 
equation of state and the computation is performed for a wide range of the polytropic index 
n (= 0.5,0.66667,0.8,1,1.25), and compactness of neutron stars {M/R)oo (= 0.03 — 0.3), where 
M and R denote the mass and radius of neutron stars in isolation. Because of the assumption of 
the irrotational velocity field, a sequence of fixed rest mass can be identified as an evolutionary 
track as a result of radiation reaction of gravitational waves. Such solution sequences are computed 
from distant detached to innermost orbits where a cusp (inner Lagrange point) appears at the inner 
edges of the stellar surface. The stability of orbital motions and the gravitational wave frequency at 
the innermost orbits are investigated. It is found that the innermost stable circular orbits (ISCO) 
appear for the case of stiff equation of state with n < 2/3. We carefully analyze the ISCO for n — 0.5 
and show that the ISCO are mainly determined by a hydrodynamic instability for {M/R)oo ^ 0.2. 
We also investigate the total angular momentum and the specific angular momentum distribution of 
the binary configuration at the innermost orbits, where the final merger process starts. From these 
quantities, we expect the final outcomes of the binary neutron star coalescence. 
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I. INTRODUCTION 

The last 10 minutes of binary neutron stars is one of 
the most important targets of laser interferometric grav- 
itational wave detectors such as LIGO , as well as a 
possible candidate for various astrophysical phenomena 
such as 7-ray bursts |^ and optical transient events [Q. 
The close binary neutron stars are in quasi circular orbits 
and evolve as a result of gravitational wave emission to- 
wards merger. The evolutionary track of binary neutron 
stars is roughly divided into three phases; (i) inspiral- 
ing phase in which the characteristic orbital separation 
a is much larger than the radius of neutron stars R and 
the binary evolves adiabatically, (ii) intermediate phase 
in which the binary still evolves in the adiabatic manner 
but a/i? is so small 2 — 5) that the effects of tidal defor- 
mation of the neutron stars become important, and (iii) 
the final coalescing phase in which the two neutron stars 
merge on the dynamical timescale. The phase (i) has 
been treated by an analytic method using a post New- 
tonian approximation and equations of motion for point 



particles, resulting in recent successful developments [||. 
On the other hand, phases (ii) and (iii) have to be treated 
numerically taking into account effects of the finite size, 
tidal deformation and hydrodynamic interaction of the 
two neutron stars. In this paper, we focus on phase (ii). 

The signal of gravitational waves in phases (ii) and (iii) 
is expected to give a wide variety of information about 
neutron stars and strong gravitational fields. In particu- 
lar, the signal around the innermost stable circular orbit 
(ISCO) at which the orbit of binaries becomes dynami- 
cally unstable to starting to merge can bring important 
information about the structure of the neutron stars. At 
the ISCO, the signal is likely to change abruptly from a 
quasi periodic wave form to a burst wave form. In the 
frequency domain, such a change will be recognized as a 
clear cliff Q . The characteristic frequency at the edge of 
the cliff may be used for extracting information about the 
nature of equations of state of high density neutron star 
matter, because the location of the ISCO may depend 
sensitively on the equations of state (or compactness) of 
the neutron stars B. The expected frequency for the 
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ISCO, /gwj can be estimated as 
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where Mq and Mq are the angular velocity, the rest 
mass of each star and the solar mass. Thus, for a typ- 
ical value of Mq ^ I.SMq, /gw may be smaller than 
IkHz, which will be in the sensitive frequency band of 
laser interferometers 0, for a relatively small value of 
MqQ. - 0.02 (see, e.g., Tables II-III). This implies that 
the characteristic frequency may be detected by laser in- 
terferometers in the near future to constrain the equation 
of state of neutron stars. One of the main purposes of 
this paper is to investigate /gw- 

Other motivations for this paper are to prepare initial 
conditions for numerical simulations of merging of binary 
neutron stars in full general relativity as well as to pre- 
dict possible outcomes in simulations starting from such 
initial data sets. Recently, several groups have been de- 
veloping numerical codes for such simulation and 
it is now becoming feasible to perform the simulation sta- 
bly and fairly accurately (see, e.g., [0). As indicated in 
[ pT[ , outcomes of the merger depend sensitively on initial 
velocity field and compactness of the neutron stars, so 
that we need to prepare the initial conditions as realis- 
tically as possible in order to obtain reliable results. It 
should also be emphasized that the simulation is not still 
an easy task and that it is very helpful to retain a reli- 
able method to cross-check the results. Careful analysis 
of the initial data sets is one of the suitable methods for 
such cross-checking. 

In this paper, we present numerical solutions for irro- 
tational binary neutron stars of equal mass in quasiequi- 
librium states. We construct such states assuming the 
existence of a helicoidal Killing vector (see Eq. (^) be- 
cause the emission timescale for gravitational waves is 
longer than the orbital period even just before merging 
and hence the system is almost in a stationary state in 
the comoving frame. We also assume an irrotational ve- 
locity field for the neutron stars because it is recognized 
as a realistic one for binary neutron stars as long as the 
spin period of the neutron stars is not as short as about 1 
millisecond [l^ . Under these assumptions, the relativis- 
tic Euler equation can be integrated to give a Bernoulli- 
type equation p^-p^ , resulting in a great simplification 
for handling the matter equations. 

We adopt a polytropic equation of state in the form 



where P, p, n and k are the pressure, rest mass den- 
sity, polytropic index, and polytropic constant, respec- 
tively. The adiabatic assumption {i.e., k is a constant) is 
appropriate because the timescale for heating and cool- 
ing inside neutron stars is much longer than the inspi- 
raling timescale of binaries due to gravitational radia- 
tion. Since realistic equations of state are not clear yet 
|l7| , we choose a wide range for n between n = 0.5 and 
1.25 so that various moderately stiff equations of state 
for neutron star matter can be approximated. We also 
survey a wide range of the compactness parameter from 
{M/R)oo = 0.03 to 0.3 where {M/R)^ denotes the com- 
pactness of neutron stars in isolation. With these param- 
eter sets, we will be able to explore the effects of various 
realistic equations of state appropriately. 

This paper is organized as follows : In Sec. ||, we 
briefly review our computational method and definition 



of physical quantities. In Sec. [II, we show numerical re- 



sults. We analyze the quasiequilibrium states from var- 
ious points of view. First, we construct a large number 
of quasiequilibrium sequences of constant rest mass for 
a wide range of n and {M/R)oo. Then, we investigate 
the stability of orbital motions for these sequences. The 
stability is analyzed by searching the simultaneous min- 
ima of energy and angular momentum as a function of 
the orbital separation (or angular velocity) along each 
sequence [^8|jl^]. The mechanism of the instability and 
the gravitational wave frequency at the innermost orbits 
are determined. We also show that the maximum den- 
sity decreases with decreasing orbital separation. Finally, 
we analyze the non-dimensional angular momentum pa- 
rameter q = Jtot/AfloM' where Jtot and Madm are the 
total angular momentum and ADM mass of the system, 
and the distribution of the specific angular momentum 
of mass elements for binaries at the innermost orbits to 



predict possible outcomes after the merger. In Sec. IV, 



we summarize the results presented in Sec. [II. In the 
Appendix, we demonstrate for completeness that neu- 
tron stars in the irrotational binary systems computed in 
this paper are stable against radial gravitational collapse. 
Throughout this paper, we take units in which G — 1 — c 
where G and c are the gravitational constant and speed 
of light. Numerical computations are performed using 
spherical polar coordinates (r, 9, ip). Greek and Latin in- 
dices run over i(or 0), r, 6, ip and r, 6, ip, respectively. 
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II. THE METHOD FOR CONSTRUCTING A 
QUASIEQUILIBRIUM SEQUENCE 

A. Outline of formulation 

Orbits of binary neutron stars shrink as a result of 
radiation reaction to gravitational waves. The ratio of 
the radiation reaction tiniescale to the orbital period is 
evaluated to be EGl 
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In deriving Eq. (j^), we adopted the quadrupole formula 
and the energy equation in Newtonian theory. Eq. (j^) 
implies that even just before the merger at a ^ 6A/adm, 
the reaction timescale is still longer than the orbital pe- 
riod and that the binary is approximately in a quasi sta- 
tionary state in the comoving frame with orbital angular 
velocity fi. Thus, neglecting the small effect of gravita- 
tional radiation reaction, we assume the existence of a 
helicoidal Killing vector in the form 



di. 



n 



d_ 



(4) 



Next, we assume that the velocity field of the neu- 
tron stars is irrotational. Viscosity inside neutron stars 
is considered to be too weak to synchronize the spin with 
the orbital rotation on the emission timescale for gravita- 
tional waves by the binary |l^ . In the final stage of inspi- 
raling, the orbital period is about 2 milliseconds. Thus, if 
the spin period of the neutron stars is much longer than 2 
milliseconds, the effect of the spin just before the merger 
is negligible so that the irrotational velocity field can 
be approximately achieved. In this assumption, further- 
more, the relativistic Euler equation can be integrated 
to give a Bernoulli- type equation [jl^p^ . Consequently, 
procedure for obtaining a solution of the hydrodynamic 
equations is considerably simplified. 

Finally, we assume that the spatial component of the 
metric is conformally flat and write the line element 
in the form |2l|-|| 

ds^ = -(a2 _ f}^i3')dt^ -h 2(3,dx'dt + -f^f.jdx'dx^, (5) 

where a, /3% \1/ and fij are the lapse function, shift vector, 
conformal factor, and flat spatial metric, respectively. 
The elliptic-type equations for a, and ^' are derived 
from the Hamiltonian constraint, momentum constraint 



and slicing condition in which Kij^^^ = with the as- 
sumption 7y = '^^fij HlHOl where Kij denotes the 
extrinsic curvature. The three metric is assumed to be 
conformally flat for simplicity. Although the solutions 
obtained in this framework are valid in general relativity 
in the sense that they satisfy the Hamiltonian and mo- 
mentum constraints, they are approximate as quasiequi- 
librium states because of the simplified formulation. We 
suspect that the solutions shown below slightly and sys- 
tematically deviate from the correct ones It should 
always be kept in mind that we need to develop a bet- 
ter theoretical framework to suppress the slight deviation 
(see [|6| for one proposal). 

Solutions of quasiequilibrium states of irrotational bi- 
nary neutron stars are computed using a numerical code 
developed by Uryii and Eriguchi [Q. The numerical 
method is described in (2^j2^ to which the reader is re- 
ferred for details. The code has been tested by comparing 
with numerical results of other groups {e.g., ||2^), as well 
as by carrying out several convergence tests [p8| . 

B. Definition of variables 

First, we define the total rest mass Mq, ADM mass 
Madm and angular momentum Jtot as 



Mo,tot= / pau^^^dV, 
Madm- 



phiauy - P + 



Jtot= J phau"u^-^'^dV, 



(6) 

^-^dy, (7) 



(8) 



where denotes the four velocity and h = l+{n+l)P/ p. 
The integral is carried out over the whole three space. 

We define the coordinate length of the semi-major axis 
Rq and half of orbital separation d as 



Ro — {Rout — ^in)/2, 
d = (i?out + Rin)/2, 



(9) 
(10) 



where Rin and i?out denote distances from the center of 
mass of the system to the inner and outer edges of the star 
respectively along the major axis. From these variables, 
we also define d = d/Ro (note that d = 1 when the 
surfaces come into contact and d oo for d oo). 
In the following, we adopt d to specify a model along a 
quasiequilibrium sequence. 
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We also define half of the separation do in another way 

as 

da = / \x\pau°^'dV, (11) 

where a; denotes a coordinate along the major axis. Here- 
after, we often refer to do as the half of the orbital sep- 
aration. 

In this paper, we compute quasiequilibrium sequences 
fixing K to be constant. This treatment is appropriate for 
the late inspiraling stage of binary neutron stars because 
the timescale of heating and cooling is much longer than 
the inspiraling timescale of the binaries due to gravita- 
tional radiation reaction. Then, the physical units enter 
the problem only through the constant n, which can be 
chosen arbitrarily or otherwise completely scaled out of 
the problem. Thus, we use non-dimensional quantities 
normalized by an appropriate power of k 



Ro 




(12) 
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= K-^/H4o,tou 


(14) 
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(15) 


Jtot 
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In the following, we often refer to a non-dimensional an- 
gular momentum parameter q = «/totA^ADM- ^'^^ 
venience, half of the total rest mass, ADM mass and 
angular momentum are also defined as Mq — Afo. tot/2, 
M = Madm/2 and J = Jtot/ 2 with their normalized 
values Mo = ^-"/^^o, M = n-'^/^M, and J = k-"J. 
Strictly speaking, M is not the mass of one star because 
the binding energy between two bodies and other inter- 
action effects are included in it. However, the fraction is 
expected to be fairly small, so that we often refer to M 
as the mass of each star. 

III. RESULTS 

A. Stability of quasiequilibrium sequences of binary 
neutron stars 

Since the vorticity is almost exactly conserved during 
the evolution of binary neutron stars, a quasiequilibrium 
sequence of irrotational binary neutron stars of constant 
rest mass and with decreasing the orbital separation (or 



increasing f?) can be identified as an evolutionary se- 
quence of binary neutron stars as a result of gravitational 
wave emission. We compute the sequences for a wide 
range of parameter space as 0.5 < n < 1.25 and typically 
as 0.1 < {M/R)oo < 0.19. (Hereafter, we use {M/R)oo 
instead of the rest mass to specify a certain sequence.) In 
some cases, we extend the analysis to {M/R)oo < 0.1 and 
{M/R)oo > 0.2. In this paper, we restrict our attention 
only to binary neutron stars of equal mass. 

Each sequence is computed gradually decreasing d 
from 3 to 1. It is found that the sequences are always ter- 
minated at an innermost orbit with d > 1 and da = da 
at which neutron stars have cusps at the inner edges of 
the stellar surfaces. This property is found irrespective 
of n and (A//i?)oo(^ 0.2). The cusps correspond to the 
inner Lagrange (LI) points. Therefore, neutron stars at 
do = da are subject to mass transfer from the cusps and 
are likely to form a dumbbell-like structure (i.e., the sys- 
tem will have a bridge between the two stars and will 
not be a "binary" any longer for da < dji). Even at 
this orbit of cusps, the timescale of gravitational radi- 
ation reaction is longer than the orbital period, except 
for extremely compact binaries, so that such a dumbbell- 
like object could be still in a quasiequilibrium state and 
may be computed with the same strategy as used here. 
However, we stop the computation at dc = da since the 
present numerical code is suitable only for binary con- 
figurations. We show several physical quantities at this 
orbit in Table I. 

To determine the dynamical stability of the orbit of 
binaries, we search for minima of J and M as functions 
of the orbital separation (or the angular velocity) along 
each sequence. Lai, Rasio, and Shapiro [Q have shown 
the existence of minima of J and M for an ellipsoidal 
model of a Newtonian irrotational binary, which indicates 
the onset of dynamical instability. (See Lombardi, Rasio 
and Shapiro [p^ for an extension including a part with 
post Newtonian corrections.) These minima of J and 
M should appear simultaneously at the same separation 
because the following relativistic identity probably holds 
for irrotational sequences with constant rest mass; 

dMADM = ^dJtot- (18) 

In our numerical computation, this identity is satisfied 
typically to < 10% except near to turning points where 
M and J are almost constant and consequently the differ- 
ence between the neighborhoods is not accurately com- 
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puted. Also, we checked that the minima of J and M 
appear simultaneously at the same separation da = c'dyn 
whenever they appear on a solution sequence ||29[] . In 
view of this, we define the location of the simultaneous 
minima of J and M as the ISCO. 

In Fig. ^ we show J as a function of da for n = 0.5, 
0.66667 and 0.8, and for (M/R)^ = 0.19. (Note that 
we have found essentially the same behavior for 0.1 < 
{M/R)ao < 0.19.) For n = 0.5, J and M have the min- 
ima at c?G = ddyn > da and so the orbits of the binary 
neutron stars will be dynamically unstable to merger be- 
fore mass transfer sets in. For n = 0.66667, the min- 
ima are located near dR and so dynamical instability and 
mass transfer will set in almost together. On the other 
hand, for n > 0.8, minima are not found for da > dj?, and 
so the mass transfer will set in before the orbit becomes 
dynamically unstable. 

Several quantities at do = ddyn for n — 0.5 and 0.66667 
are summarized in Tables II and III, respectively. As 
shown in Table I, d at dc = da depends very weakly 
on {M/R)oo- In contrast, d dX da = ddyn becomes 
larger with increasing {M/R)oo, as shown in Tables II 
and III. For n = 0.5, an ISCO always exists for any 
{M/R)oo < 0.3. For n = 0.66667, no ISCO exists for 
the sequence of {M/R)oo < 0.17, but there is one for 
{M/R)oo > 0.17. (Note that the curves of J and M 
as functions of da suggest that da is close to da even 
for {M/R)oo < 0.17 for n = 0.66667.) The results for 
the smaller compactness {M/R)oo < 0.17 with various 
n are consistent with the result in the Newtonian limit 
{M/R)oo — > |2^, implying that the above property of 
the dynamical instability holds irrespective of compact- 
ness as long as {M/R)oo < 0.17. 

To summarize, we show a schematic figure in Fig. |^. 
We have found two cases for evolution of close binary neu- 
tron stars with irrotational velocity field. In one case, the 
sequence is terminated when cusps appear before the bi- 
naries reach the minima of J and M as shown in Fig. H(a). 
In this case, mass transfer will occur before the orbit be- 
comes dynamically unstable. Binary neutron stars with 
n > 2/3 are in this category. In the other case, the bina- 
ries reach the minima before the cusps appear as shown 
in Fig. 11(b). In this case, the orbit becomes dynamically 
unstable before mass transfer sets in. Binary neutron 
stars with n < 2/3 are in this category. Based on this 
result, we hereafter refer to the orbit at da = ddyn for 
n ^ 2/3 and at da = dji for n > 0.8 as "innermost orbit" . 



B. Frequencies of gravitational waves at the final 
orbit 



As discussed in Sec. I, it is important to clarify the fre- 
quency of gravitational waves at the ISCO from the view- 
point of gravitational- wave-astronomy. For n = 0.5 and 
0.66667, we find that an ISCO (i.e., a point at which dy- 
namical instability sets in) does exist, and the frequency 
is computed from Tables II and III as 



/gw 



0.8 



'1.5Mq 



1.1 - 1.15 



1.35-1.4 



Mo 
'1 .5Mq ' 

Mo 



for (M/i?) oo = 0.14, 
[kHz for (M/i?)oo = 0.17, (19) 
I kHz for {M/R)oo = 0.19. 



We note that the ISCO here for {M/R)^ < 0.2 is prob- 
ably determined by the hydrodynamic instability JTsf , 
but not by the general relativistic orbital instability. The 
reason is that the general relativistic orbital instability 
should not depend on {M/R)oo in the absence of hy- 
drodynamic effects |3^. Even in the presence of such 
effects, A/adm^^ should depend only weakly on {M/R)oo 
when the general relativistic orbital instability dominates 
to determine the ISCO. However, for (M/R)^ < 0.2, 
Mabm^ at the ISCO does depend strongly and system- 
atically on {M/R)oo (cf, Tables II and III, and Fig. |). 



We will discuss this point in detail in Sec. HI C . 

For a realistic neutron star with Mq ~ 1.5A/0, the ra- 
dius will be in the range between 10 and 15km. This 
implies that (M/R)^ is between 0.14 < {M/R)oo < 0.2 
for a typical ADM mass ~ 1.4AfQ where Moo de- 
notes the ADM mass in isolation. The present result 
suggests that if the neutron star radius is fairly large, 
i.e., ~ 15km, the frequency of gravitational waves at the 
ISCO will be less than IkHz which is in the sensitive 
frequency band of the laser interferometers. However, if 
the radius of neutron stars is 10 km, the frequency 
at the ISCO is larger than IkHz and it will be difficult 
to detect gravitational waves emitted from binary neu- 
tron stars with Moo ~ 1.4M0 orbiting near the ISCO by 
means of the interferometers [Q. 

For n — 0.8 and 1, the frequency of gravitational waves 
at da = dfi is given approximately by Eq. (19), implying 
that /gw would be larger than the frequency shown in 
Eq. (|l|) for these cases. Hence, even for R ~ 15km and 
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Mq ~ I.5M0, the frequency at ISCOs may be larger than 
IkHz. 

It should be pointed out that for n > 0.8, mass trans- 
fer is likely to set in before the binaries reach the ISCO 
as discussed above. If the timescale of mass transfer 
is longer than the timescale for evolution due to grav- 
itational wave emission, the signal of the gravitational 
waves may not change abruptly from a quasi periodic 
type to a different one at this orbit. However, if the mass 
transfer proceeds quickly, a characteristic signal may be 
produced. As mentioned above, the frequency of gravi- 
tational waves at the onset of mass transfer is approxi- 
mately given by Eq. ([l^) for n = 0.8 and 1. Thus, for the 
case of relatively soft equations of state, the character- 
istic signal at mass transfer may be detected for binary 
neutron stars of {M/R)oo ^ 0.14 and Moo ~ IAMq. 

C. Origin of the ISCO 

In Fig. ^, we plot Madm^^ at ISCOs for n = 0.5 as 
a function of {M/R)oo using data tabulated in Table II. 
As shown in Fig. ^ and as argued in Sec. IIIB, Madm^^ 
at the ISCOs increases with increasing {M/R)oo- This 
strongly indicates that the ISCO is determined by the 
hydrodynamic instability p^ ]. On the other hand, in 
the limit {M/R)ao 1/2, hydrodynamic effects should 
become less important so that the general relativistic or- 
bital instability can determine the ISCO. In our present 
formulation, not only hydrodynamic but also general rel- 
ativistic effects (in the framework of the conformal flat- 
ness approximation) are taken into account, which en- 
ables us to investigate the origin of the ISCO. In this 
subsection, we carry out a detailed analysis using the 
results of n = 0.5. 

As a first step, we derive fitting formulae for {M/R)oo 
< 0.2 to clarify the behavior of MqH. and Mfl at IS- 
COs as functions of {M/R)oo- In the derivation, the 
following two points are taken into account: (1) In the 
Newtonian limit (M/i?)oo 0, Mo^l{M/R)^^^ and 
Mri(M/ R)o^^'^ converge to an identical constant because 
they are non-dimensional and should be identical in New- 
tonian gravity. (2) In the Newtonian limit, the ISCO is 
determined by the hydrodynamic effect jl^, and this is 
the case for the finite {M/R)oo as long as {M/R)oo is 
small. Consequently, all of the general relativistic cor- 
rections should appear in the power series of {M/R)oo 



from the viewpoint of the post Newtonian approxima- 
tion. From (1) and (2), we fix the functions of and 
Affi in the form 



Mo^{MlR)^'^ = a + b{M/R)oo + c{M/R)l^, (20) 
Mr!(M/i?)-3/2 = a' + b'{M/R)oo + c'{M/R)l^, (21) 



(namely, we expand them up to second post Newto- 
nian order) and determine the coefficients (a, 6, c) and 
(a', 6', c') by the least square fitting. We here note that 
although we know the value for a — a' from the New- 
tonian computation [2^ , we do not use the data for the 
fitting. Instead, we compute the value at {M/R)oo = 
from the fitting formula and compare with the Newtonian 
results for a cross-checking. 

In the least squares fitting, we use the data sets for n — 
0.5 of 0.03 < {M/R)oo < 0.17 (or 0.19) shown in Tables 
II. The resulting coefficients for the fitting formulae are 
tabulated in Table IV. We may roughly conclude. 



a = 0.286 ±0.001, 6 = 0.40 ± 0.03, 
a' = 0.286 ±0.001, 6' = 0.22 ± 0.02. 



(22) 



For c and c', on the other hand, a large uncertainty of 
0(0.1) exists. In a Newtonian analysis, MoQ{M/R)^-^ 
= Mn{M/R)^-^ ~ 0.284 for n = 0.5 13, which agrees 
well with a and a' in the above fitting formulae. This 
shows the validity of the fitting method adopted here. 

In Fig. |, Mon{M/R)^-^ and Mat,m^{M / R)^-^ at 
the ISCO are plotted as a function of {M/R)oc together 
with results from the fitting formulae derived above. The 
filled circles and triangles denote the numerical results, 
and the solid and dotted lines denote results obtained 
from the fitting formulae. These figures indicate that the 
fitting formulae agree fairly well with numerical results 
for (M/i?) 00 ^ 0.2, implying that the post Newtonian 
expansion is reasonable in this regime. In particular, the 
fitting is very good around {M/R)oo ^ 0.15. In this 
regime, Moil and A/adm^I are well approximated by only 
the first post Newtonian correction, i.e., the (M/R)'^ 
term is not very important. 

The results of the fitting imply that the frequency 
of the ISCO depends strongly and systematically on 
{M/R)oo- Taking into account that the ISCO in the 
Newtonian limit is determined by the hydrodynamic in- 
stability, we may naturally infer that the ISCO for small 
{M/R)oo ^ 0.2 is still determined by the hydrodynamic 
instability with general relativistic corrections. As we 
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mentioned in Sec. [II B, the compactness of a realis- 



tic neutron star is likely to be in the range between 
0.14 < {M/R)oo < 0.2. Therefore, the ISCO of the 
realistic binary neutron stars is determined by the hy- 
drodynamic instability but not by the general relativistic 
instability. 

For 0.2 < (M/i?)oo < 0.3, on the other hand, the 
fitting formulae do not agree very well with the numerical 
results, implying that the angular velocity at the ISCOs 
seems to be affected by strong nonlinear effects of general 
relativity. In this regime, the hydrodynamic effects seem 
to be less important. In Fig. ^, we show the ratio of semi- 
diameters for neutron stars at the ISCOs as a function 
{M / R)oo for n = 0.5. The semi-diameters (ai, 02, 03) are 
computed from 



(23) 



a„ = / ^^di (n = 1,2,3), 



where the line integral is taken along straight lines in 
Cartesian coordinates, and the three directions (n = 
1 — 3) are chosen to be orthogonal each other, oi is the 
longest diameter along a line connecting the centers of 
the two neutron stars. 02 and 03 are those in the equa- 
torial and the meridional plane, respectively, and they 
intersect each other at the coordinate center of ai . The 
circles and crosses in Fig. || denote 02/01 and a^/ai, re- 
spectively. These ratios appear to converge towards unity 
with increasing {M/R)oo- This indicates that the tidal 
effect is less important for determination of the ISCO for 
{M/R)oo ~ 0.3. However, even in this highly general rel- 
ativistic regime, Madm^^ still does not converge to a con- 
stant as shown in Fig. ^. Since Madm^^ depends strongly 
on {M/R)oo, the effect of the hydrodynamic instability 
still appears to dominate over the general relativistic or- 
bital instability for determining the ISCO. 

Cook 1^ and recently Baumgarte have investi- 
gated the ISCO for binary black holes using the confor- 
mal flatness approximation for the geometry as we do. 
According to their results, Mfl at the ISCO is 0.17- 
0.18. {A4 is the sum of the gravitational masses of the 
two black holes, and is slightly different from A/adm al- 
though the difference is not important here.) As shown in 
Fig. |, Madm^ ^ 0.11 at the ISCO for (M/R)^ = 0.3 
and it is still monotonically increasing. This tendency 
does not contradict their results. The ISCO in their 
work is purely determined by the general relativistic or- 
bital instability. Therefore, we may expect that the 



ISCO of binary neutron stars is eventually determined 
by the general relativistic orbital instability for a suffi- 
ciently large {M/R)oc > 0.3 and A/adm^^ converges to a 
constant ^ 0.17 — 0.18, although such a large compact- 
ness is unrealistic for neutron stars. (We have the value 
Madm^ — 0.18 for n = 0.5 case, if we naively extrapo- 
late data points at the ISCO with {M/R)^ = 0.24 and 
0.3 linearly to {M/R)^ = 4/9.) 

Finally, we compare the present results with ones com- 
puted using other analytic methods. Kidder, Will and 
Wiseman Damour, Iyer and Sathyaprakash |Q, and 
Buonanno and Damour ||3^ have determined the ISCO 
using relativistic equations of motion for two point parti- 
cles in which some general relativistic corrections have 
been taken into account. They have derived the or- 
bital angular velocity at the ISCO as M.VI ~ 0.0605, 
0.08850 and 0.07340, respectively. Since equations of 
motion for point particles have been used in their meth- 
ods, the ISCO has been determined by the general rel- 
ativistic orbital instability and AA^l has been indepen- 
dent of {M/R)oo- This is in contrast with our results, 
in which Madm^^ depends on {M/R)ao even around 
{M/R)oo ^ 0.3 and consequently the general relativis- 
tic orbital instability does not seem to dominate. 

The reason for this contradiction is clear for {M/R) 00 
< 0.2. As we showed above, the ISCO for iM/R)^^ < 0.2 
is determined by the hydrodynamic instability. In this 
regime for the compactness, Madm^^ ^5 0.06 which is 
smaller than their results, namely, that before the gen- 
eral relativistic effects between two stars become impor- 
tant, the orbits become unstable because of the hydro- 
dynamic instability. Since realistic neutron stars have 
0.14 < {M/R)oo ^ 0.2, the approaches using equations 
of motion for point particles are not appropriate for in- 
vestigating the ISCO of binary neutron stars. 

For 0.2 < {M/R)oo < 0.3 and Madm^^ < 0.11, 
our present results suggest that the ISCO may be de- 
termined by the hydrodynamic instability. According 
to the point particle approaches, however, the general 
relativistic orbital instability should be important for 
Maum^ > 0.06 — 0.09 implying that our results contra- 
dict theirs in this regime, too. The reason for this con- 
tradiction is not clear. The works of ourselves. Cook and 
Baumgarte are based on the conformal flatness approxi- 
mation which may fail to include some important general 



relativistic effects 25 1. The neglected effects might play 
an important role in determining the ISCO for the regime 
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{M/R)oo > 0.2. Thus, we should not draw any strong 
conclusions for such highly relativistic binaries from re- 
sults obtained in this framework since this might be the 
reason for the contradiction. To pin down the uncertainty 
for the high compactness regime, we need to carry out 
computations including higher general relativistic correc- 
tions. 

D. Evolution of the maximum density 

In Figs. H(a)-(c), we show the relative change in the 
maximum energy density of a star Ae = (cmax — eoo)/eoo 
as a function of d for n = 0.5,0.66667 and 0.8, respec- 
tively. Here, e ~ p + nP, and Cmax and Coo denote re- 
spectively its maximum value for d < oo and its value for 
an isolated spherical star {d = oo) as computed from the 
TOV equation. 

The figures indicate that the maximum density slightly 
decreases with decrease of the orbital separation irrespec- 
tive of n for c? < 2. For larger separations of d > 2, each 
star becomes almost spherical so that Ae can be less 
than the numerical uncertainty. The small global shifts 
(< 1%) of each curve from zero for 2 < d < 3 are thought 
to be due to numerical errors associated with our numer- 
ical scheme (a finite difference scheme). Small deviations 
around the average for 2 < d < 3 are due to a system- 
atic error associated with the Legendre expansion for the 
gravitational field. (See for details.) Thus, our com- 
putation indicates no evidence that neutron stars in bi- 
nary systems become dynamically unstable against grav- 
itational collapse. Instead, the maximum density slightly 
decreases with decreasing orbital separation irrespective 
of n. (See the Appendix for more detailed computations 
for selected values oi n.) This result is consistent with 
recent analytic calculations in general relativity |36| . 

E. Prospects for the outcome after the merger 

First, we investigate the quantity q = Jtot/Mj^-^-^ 
which is an important parameter for predicting the out- 
come after a merger. If the q of the object resulting from 
the merger is larger than unity, it will not be able to col- 
lapse to form a black hole. Thus, if q is larger than unity 
for a binary at its innermost orbit, there is a possibil- 
ity that no black hole will be formed promptly after the 
merger. On the other hand, if q is less than unity before 



the merger, we could consider that the system is a candi- 
date to become a black hole, because q will not increase 
during the merger due to emission of gravitational waves 

In Fig. ^, we show the relation between {M/R)oc and 
q at the innermost orbits of irrotational binary neutron 
stars for n — 0.5 — 1.25 and we also show the relation for 
corotating binary neutron stars with n = 1 for compari- 
son. Here, the innermost orbits for irrotational binaries 
are the orbits at da = ddyn for n — Q.b and at do = djf for 
n > 2/3 as we define in Sec. HI A . (Note that Hr ~ ddyn 
for n = 2/3.) For corotating binaries, we choose the or- 
bit at the energy and angular momentum minima. It is 
found that q is always less than unity irrespective of n 
for irrotational binaries with {M/R)oo > 0.13. As men- 
tioned in Sec. Ill A, {M/R)oo will be larger than 0.14 
for a realistic neutron star. Therefore, q for realistic bi- 
nary neutron stars decreases below unity before merger, 
implying that the binary system satisfies a necessary con- 
dition to form a black hole. 

Here, two cautions are appropriate. (1) For very stiff 
equations of state, the total rest mass of the system can 
be smaller than the maximum allowed rest mass of a 
spherical star. Indeed, for n = 0.5 and {M/R)oo = 0.14, 
this is the case. Thus, for very stiff equations of state, a 
black hole will not be formed even if g < 1 before merger 
and {M/R)oc is fairly large (~ 0.14). (2) Even in the case 
when the system is supramassive in the sense that the rest 
mass of the system is larger than the maximum rest mass 
of a spherical star with the same equation of state, the 
condition q < 1 is not sufficient for black hole formation, 
because rapid rotation could support the self-gravity of 
the supramassive star [ p8[ . To find the outcome of the 
merger, numerical simulations are obviously necessary. 

In contrast with the above results, q is larger than unity 
even for {M/R)oo ~ 0.17 at the energy and angular mo- 
mentum minima for corotating binary neutron stars with 
n — 1 The reason for the difference between the 

two cases is that stars in the irrotational binaries have 
only negligible spin angular momentum while those in the 
corotating binaries have significant spin angular momen- 
tum. In the study of corotating binary neutron stars with 
{M/R)oo ^ 0.17, one could reach the conclusion that a 
black hole might not be formed after the merger, which 
is completely different from our present conclusion for ir- 
rotational binary neutron stars. This result brings us to 
another caution, that if we adopted corotating velocity 
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fields, we could reach an incorrect conclusion for realistic 
binary neutron stars which will have (nearly) irrotational 
velocity fields. 

Next, we investigate the mass spectrum with respect 
to the specific angular momentum of fluid elements in 
order to estimate the mass of the disk around a black 
hole formed after merger. We consider this for binaries 
at innermost orbits which are chosen in the same manner 
as above. The mass spectrum with respect to the specific 
angular momentum is defined as 

Mo Mo Jj,yj 

where the integral is performed for fluid elements which 
have specific angular momentum larger than a given 
value j. The specific angular momentum j is defined 
as 

j = hu^. (25) 

Here, we note that j of each fluid element is exactly con- 
served in axisymmetric systems for ideal fluid. For non- 
axisymmetric systems, it may not be conserved due to 
some dissipation processes such as gravitational radia- 
tion or other outward transportation process. 

In Figs. H, we show Mo(j)/Mo as a function of j/MAOM 
for n = 0.5 - 1.25, and (M/R)^ = 0.14 and 0.17. Here, 
we choose binaries with {M/R)oo > 0.14 in order to 
be realistic. The vertical dotted lines in Figs. || denote 
j/AfADM of a test particle orbiting at the ISCO around a 
Kerr black hole for q = 0.99, 0.95, and 0.9 showing 
the minimum allowed values for test particles orbiting 
a black hole. For n = 1 (Fig. ||(d)), we also show the 
results for corotating binaries for comparison. The fig- 
ures show that the mass fraction of fluid elements with 
j'/Madm > 1-6 for irrotational binary neutron stars is 
zero irrespective of n and {M/ R)oo- In contrast, that for 
corotating binaries is ~ 5% even for j/A/adm — 2. 

Now, we assume the following hypotheses which are 
likely to hold in the merger of binary neutron stars with 
{M/R)oo > 0.14: (1) A black hole is formed. This as- 
sumption could be wrong for n — 0.5 and {M/R)oo ~ 
0.14, but would be correct for a moderate stiffness, 
2/3 ^ n 1 (for n = 1, numerical simulations have 
indicated the correctness of this assumption |jri|). (2) 
The disk mass formed after the merger is much smaller 
than the black hole mass (say at most 10% of the total 
rest mass), and consequently the mass of black hole Mbh 



is nearly equal to the initial gravitational mass Madm 
(say > 0.9Madm)- (3) The specific angular momentum 
j for most of fluid elements changes by at most 10% dur- 
ing the merger process. One mechanism to change j is 
gravitational radiation which would decrease j of all fluid 
elements by < 10%. The other is hydrodynamic interac- 
tion which could transport j outward and consequently a 
part of fluid elements could have larger j than its original 
value. Therefore this hypothesis would be reasonable if 
the merger process to form the black hole proceeds on the 
dynamical timescale of the system and it is shorter than 
the timescale for the outward transportation of angular 
momentum by the hydrodynamic interaction ]4C| ]. (4) 
The quantity q of the black hole finally formed is about 
10% smaller than the initial value for the system because 
of gravitational radiation reaction {i.e., q ^ 0.9). 

With the hypotheses (1), (2) and (4), the specific an- 
gular momentum of a test particle at the ISCO around 
the formed black hole is approximately determined from 
the Madm and Jtot given initially. As shown above, q for 
binaries with {M/R)oo > 0.14 is less than unity and ac- 
cording to (4), q of the black hole will be ^ 0.9 after the 
merger. This implies that the specific angular momen- 
tum j at the ISCO around the black hole would be larger 
than 2Mbh |4^, and to form disks, there have to exist 
a fraction of fluid elements of j > 2Mbh ^ 1-8Madm- 
However, as shown in Figs. |[ the mass of the fluid ele- 
ment of j/A/adm > 1-6 is almost zero before merger for 
any n. Thus, if all the hypotheses are correct, the mass 
fraction of the disk around the black hole is almost zero 
after the merger of irrotational binary neutron stars of 
{M/R)oo > 0.14. In contrast, the mass fraction of the 
disk with j > 2Mbh can be about ~ 5% for corotating 
case as shown in Fig. ||(d). 

Recent numerical simulations of merging binary neu- 
tron stars of equal mass in full general relativity for n ~ 1 
|]lT| indicate that for irrotational binaries, the mass of the 
disk around the black hole formed after merger is negligi- 
ble, while for corotating binaries, it could be '--^ 5 — 10%. 
These results are in complete agreement with the above 
analysis. This suggests that the above working hypothe- 
ses are approximately correct. Also, it suggests that the 
present analysis can be applied for other n. Thus, ir- 
respective of n and {M/R)oq{> 0.14), a massive disk is 
unlikely to be formed around the black hole after mergers 
of irrotational binary neutron stars. 

The present result leads us to give a strong caution 
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with regard to previous numerical studies. As shown 
above, the fraction of the disk mass depends sensitively 
on the initial velocity field and effect of general relativity 
(which determines the minimum allowed value of the spe- 
cific angular momentum of the disk) . In the past decade, 
most of simulations of binary neutron stars of equal mass 
have been performed using rather crude initial conditions 
for the velocity field and without general relativistic ef- 
fects (see, e.g., for review). Such simulations in some 
cases have provided results in which a disk is formed 
around the merged objects. However, the present analy- 
sis suggests that such results might be wrong and gives 
a warning that it is dangerous to draw any conclusion 
about the mass of the disk without fully including gen- 
eral relativistic effects and without using correct velocity 
fields for the initial conditions. 

IV. SUMMARY 

We have presented numerical results for quasiequilib- 
rium sequences of irrotational binary neutron stars of 
equal mass in general relativity adopting the polytropic 
equation of state p — Computations have been 

performed for a wide range of the polytropic index n 
and the compactness {M/R)oo. We have analyzed the 
sequences from various points of view. The conclusions 
obtained in this paper are as follows. 

The dynamical instability for orbits of irrotational bi- 
nary neutron stars sets in before mass transfer occurs for 
n < 2/3 (r > 2.5). Thus, the binary reaches the ISCO 
without mass transfer in this case. On the other hand, for 
n > 2/3 (r < 2.5), mass transfer will occur before bina- 
ries reach the ISCO. This property depends very weakly 
on the compactness for {M/R)oq < 0.17. The frequency 
of gravitational waves at the innermost orbit (ISCO or 
the point where mass transfer sets in) is about 800 to 1500 
Hz depending on the compactness. For a binary of less 
massive {i.e., less compact) neutron stars, the frequency 
could be less than 1000 Hz and as a result the charac- 
teristic signal of gravitational waves emitted around the 
ISCO may be detected by laser interferometric detectors. 

We have determined the ISCOs for a wide range of 
{M/R)oo < 0.3 for n = 0.5 and of 0.17 < (M/i?)oo < 
0.27 for n = 0.66667. The results indicate that the ISCO 
is determined by the hydrodynamic instability and not 
by the general relativistic orbital instability for realistic 



binary neutron stars with 0.14 < {M/R)oo ^ 0.2. Our 
present results indicate that the approaches using equa- 
tions of motion for point particles are not appropriate for 
determining the ISCO of binary neutron stars if the neu- 
tron stars are not extremely compact {M/R)oo 3> 0.2. 

The maximum density of neutron stars in irrotational 
binary systems decreases as the orbital separation de- 
creases. Associated with this, we show that neutron stars 
in irrotational binary systems (which are stable in isola- 
tion) are stable against gravitational collapse until the 
merger sets in (see Appendix for details). 

The angular momentum parameter q — ^tot/AfloM 
is less than unity for {M/R)oo ^ 0.13 irrespective of 
n for irrotational binaries at the innermost orbits {i.e. 
da = or ddyn)- Since {M/R)oo of realistic neutron 
stars is expected to be larger than 0.14, the result im- 
plies that realistic binary neutron stars satisfy a neces- 
sary condition for forming black hole, q < 1, before the 
merger. However, this condition is not sufficient for de- 
termining the outcome after the merger as mentioned in 



Sec. pi E\ Numerical simulations in full general relativ- 
ity are obviously necessary to clarify the nature of the 
object resulting from the merger pT[ |. 

It is found that the specific angular momentum of all 
of the fluid elements in irrotational binary neutron stars 
of equal mass will become very small (j'/Madm < 1-6) 
before merger irrespective of n and {M/R)oo- As long as 
outward transportation of angular momentum does not 
work effectively during merging, the disk mass around a 
black hole which would be formed after the merger for 
{M/R)oo ^ 0.14 will be negligible. We should note here 
that the present conclusion is obtained for binaries of 
equal mass (or nearly equal mass). In the case when the 
mass ratio of two stars deviates from unity by a large fac- 
tor, the star of smaller mass could be tidally disrupted 
by the companion of larger mass at a fairly large orbital 
separation at which the specific angular momentum of 
fluid elements could be large enough to form a disk dur- 
ing the tidal disruption. For exploring the possibility for 
formation of a disk in this way, it is necessary to perform 
computations for binary neutron stars of unequal mass, 
and this is one of the key issues for the future work. 
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APPENDIX: STABILITY OF NEUTRON STARS 
AGAINST GRAVITATIONAL COLLAPSE AND 

EXISTENCE OF SUPRAMASSIVE 
IRROTATIONAL BINARY NEUTRON STARS 

It is important to check that neutron stars which are 
stable in isolation are stable against gravitational col- 
lapse to black holes even in binary systems. (See [ pT|j23| 
about a claim that the neutron stars in irrotational bi- 
nary systems may be unstable to collapse individually.) 
In this Appendix, we show that no evidence is found in 
our numerical experiment for selected values of n. 

The concern about gravitational collapse is relevant 
for neutron stars of near to the maximum allowed rest 
mass. Thus we here focus only on such massive neutron 
stars. As noted in the caption of Table I, the compactness 
{M/ R)oo of neutron stars at the maximum allowed rest 
mass is larger than 0.25 for n < 0.8, ^ 0.21 for n — 1 and 
~ 0.17 for n — 1.25. The compactness of realistic neutron 
stars would be in a range between 0.14 < {M/R)ao ^ 
0.2 and probably less than 0.25. From this fact, we pay 
attention only to n = 1 and 1.25 here. 

First, for calibrating our numerical code, we compare 
numerical results at large d with those for spherical stars 
computed by solving the TOV equations. Then, we show 
the existence of supramassive neutron stars in close bi- 
nary systems confirming the excess of the maximum rest 
mass beyond the maximum rest mass of the spherical 
stars. Finally, we discuss the stability against collapse of 
supramassive neutron stars in close binary systems (stars 
whose rest mass is larger than the spherical maximum 
mass) . 



In the analysis, we computed quasiequilibrium states 
with fixed d, changing {M/R)^ (or pmax)- In Fig. |, 
we show Mq as a function of Pmax (solid lines) (a) for 
n ~ 1 and (b) for n = 1.25. Here, d is chosen between 
1.3125 and 2.0 for n = 1 and between 1.3125 and 2.5 
for n = 1.25. For obtaining the configuration at d = oo 
(denoted by the dashed lines), the TOV equations are 
solved. Note that the maximum mass for spherical stars 
is Mo ^ 0.180 at Pmax ^ 0.318 {{M/R)oo ^ 0.214) for 
n = 1 and Mo ^ 0.2167 at p,nax ^ 0.148 {{M/R)^ ~ 
0.172) for n=1.25. 

For large separation d > 2.0, neutron stars in irrota- 
tional binary systems become almost spherical irrespec- 
tive of n. Accordingly, the solid curves should converge 
toward the dashed curve with increasing d. Figure ^ 
shows that this is approximately the case, but the curve 
in the limit d ^ 1.0 slightly deviates from the dashed 
curve. This deviation is due to a systematic error in 
numerical computation. This shows that the present nu- 
merical code overestimates the rest mass by 0.4% for 
a quasiequilibrium state with a given Pmax- For this rea- 
son, the curve in the limit c? ^ 1.0 denotes the relation 
between the rest mass and Pmax for spherical stars com- 
puted in our numerical code, and hereafter, we refer to 
this curve as the relation for spherical stars. The max- 
imum of this curve is Mq™^'^ = 0.1802 for n — 1 and 
= 0.2174 for n 1.25. 

With decreasing d, the maximum value of Afo(pniax) 
increases. This suggests that the maximum allowed rest 
mass of neutron stars in irrotational binary systems is 
increased due to the tidal effect. The maximum at 
d ^ 1.3125 is Mq ~ 0.1808 for n ^ 1 and 0.2179 for 
n = 1.25. This implies that supramassive neutron stars 
exist in irrotational binary systems, but the maximum 
mass Mq can be increased by at most ~ 0.3% even at 
da — dn. This increase is much smaller than that in the 
corotating binaries ||3^, for which Mo increases by ~ 2% 
for n = I and and by ^ 1.5% for n = 1.5. The reason 
is that neutron stars in irrotational binary systems have 
negligible spin, while those in corotating binary systems 
have significant spin and the centrifugal force by this can 
increase the maximum allowed rest mass. Such binary 
systems with supramassive neutron stars could be formed 
in principle after a supramassive core collapse to undergo 
bifurcation. 

Baumgarte et al. |^ have discussed the stability of 
a supramassive neutron star in a corotating binary us- 
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ing the turning point method . They have concluded 
that the supramassive and aU other neutron stars below 
a critical density in a corotating binary are stable against 
gravitational collapse. Following their discussion, we plot 
in Fig. |o| sequences of constant Mq for a variety of Mq 
on the Pma^-M plane (the solid and thick long dashed 
lines). Each solid curve is drawn for 1.3125 < d < 2.0 in 
the case n = 1 and for 1.3125 < d < 1.875 in the case 
n — 1.25. Here, each sequence of constant Mq can be 
interpreted as an evolutionary sequence of a binary sys- 



tem as discussed in section [IF When a binary system 
evolves as a result of gravitational radiation, the gravita- 
tional mass decreases. Thus, if Pmax decreases (increases) 
with decreasing gravitational mass, the neutron stars are 
expected to be in a stable (unstable) branch. Therefore 
we may conclude that the curves in region (i) are stable 
while those in region (iii) are unstable. 

For the thick dashed curves in region (ii), a turning 
point exists at a maximum M. In contrast to the coro- 
tating binary case, it is not fully established whether the 
stability of a solution sequence changes at this turning 
point or not. However, from the tendency of solution 
the sequences (i), (ii) and (iii), it is reasonable to expect 
that a part of the solutions on the left branches of the 
thick dashed curves in region (ii) are stable against grav- 
itational collapse. From these results, it is reasonable to 
conclude that stable supramassive neutron stars can exist 
in irrotational binary systems as a result of tidal defor- 
mation, although the excess of the rest mass beyond the 
maximum rest mass of the spherical stars is tiny. 
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TABLE I. Characteristic quantities for irrotational binary neutron stars for configurations at da = dn. Here, pmax is the 
maximum rest mass density. Note that for spherical stars, the maximum value of the normahzed rest mass and the compactness 
at the maximum (Mo, (M/R)^) are (0.152, 0.316) for n = 0.5, (0.155, 0.278) for n=0.66667, (0.162, 0.252) for n = 0.8, (0.180, 
0.214) for n=l, and (0.217, 0.172) for n=1.25. Because of an uncertainty for determination of the configuration with cusps, 
MoO and Madm^^ have an error of 1 — 2% for each (M/R)ao. 
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TABLE II. The same as Table I but for a configuration with n = 0.5 at da = ddyn- Wo include Madm^^ in this Table. 
Because of an uncertainty for determining the minima of Jtot and Madm, Mq^I and Madm^^ have an error of 1 — 2% for each 
{M/R)oo. 
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TABLE III. The same as Table II but for a configuration with n = 0.66667 at da = ddyn. Because of an uncertainty for 
determining the minima of Jtot and Madm, MqO, has an error of 1 — 2% for each {M/R)oo- Note that ddyn is nearly equal to 
dR for (M/R)^ < 0.19 but not for (M/R)^ > 0.2. 
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(0.2864,0.200,0.18) 



TABLE IV. Coefficients of the fitting formula Eq. (|2C 
for the results for n = 0.5 (Table II). The column headed 
{M/R)oo shows the data sets used for each fitting. 
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FIG. 1. J as a function of da for (a) n — 0.5, (b) 
n = 0.66667 and (c) n = 0.8 and for {M/R)oo = 0.19. The 
points of smallest da along the curves correspond to dc ~ dn 
where neutron stars have cusps at the inner edges of the stellar 
surfaces. 
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(a) n>n„ (b) n<n, 




n^j=:2/3 



FIG. 2. Schematic figure of J and M as functions of da 
for sequences of constant rest mass, dn corresponds to the 

point where neutron stars have cusps at the inner edges of the 
stellar surfaces, ddyn corresponds to a turning point along the 
sequence. 
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FIG. 4. Mo SI (filled circles) and MQ (filled triangles) at 
the ISCO {dc — rfdyn) as a function of (M/R)oo for n = 0.5. 
The solid and dotted lines denote the fitting formulae derived 
using the data sets with 0.03 < (M/R)^ < 0.17 (sohd hnes) 
and 0.03 < {M/R)ao < 0.19 (dotted lines), respectively. The 
squares at {M/R)ac, = 0.0 denote the results from the New- 
tonian computation. 




0.10 



0.15 



0.20 
(M/R)„ 



0.25 



0.30 



FIG. 3. Plot of (M/R)oo SIMadm relations at the 
ISCO {da ~ ddyn) for binary neutron stars with n = 0.5 
and n = 0.66667. 
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FIG. 5. The ratios 02/01 (circles) and 03/01 (crosses) 
as functions of {M/R)ao at ISGOs for n = 0.5. Values at 
{M/R)oo = 0.0 are those from the Newtonian computation. 
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FIG. 7. Relation between (M/R)^ and •/tot./Af^Qi^i at 
the innermost orbits of binary neutron stars. The solid lines 
denote the relation for irrotational binaries with n = 0.5, 
0.66667, 0.8 and 1 from top to bottom, and the dashed line is 
for corotating binaries with n = 1. For the irrotational case, 
the models at da ~ da arc chosen as the innermost orbits, 
and for the corotating case, the models at the energy and 
angular momentum minima are chosen. 



FIG. 6. ( )/eoo as a function of d for (a) n = 0.5, 

(b) n = 0.66667 and (c) n = 0.8. Cusps appear at the inner 
edges of the stellar surfaces at the point of smallest separation 
on each curve. 
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FIG. 8. Fraction of rest mass with specific angu- 
lar momentum larger than j for the innermost orbits for 
{M/R)^ = 0.14 and 0.17 and for (a) n = 0.5, (b) 
n = 0.66667, (c) n = 0.8, (d) n = 1, and (e) n = 1.25. 
The solid and dashed lines correspond to the results for ir- 
rotational and corotating binaries, respectively. The lines for 
{M/R)oo ~ 0.17 appear as upper linos in each panel. For 
the irrotational case, the models at da ~ ddyn are chosen 
for n — 0.5, and the models at da = da are chosen for the 
other n. For the corotating case, the models for which the 
two surfaces come into contact are chosen. 
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FIG. 9. Mo as a function of pmax (a) Solid lines from top 
to bottom correspond to d = 1.3125, 1.375, 1.5, 1.625, 1,75, 
1,875 and 2 for the n = 1 case, (b) Solid lines from top to 
bottom correspond to d = 1.3125, 1.5, 1.75,2 and 2.5 for the 
n = 1.25 case. A dashed line corresponds to a solution of the 
TOV equations in each panel. 
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FIG. 10. Sequences for constant Mo arc plotted in the 
Pmax — M plane (thick solid lines and thick long dashed lines) . 
(a) n = 1 case. Prom left to right, the Mo of each curve 
changes from 0.1797 to 0.1802 in steps of 0.0001 (solid lines) 
in region (i) (i.e., pmax < 0.31). In region (ii), two long dashed 
lines with Mq = 0.1803 and 0.1804 (upper one and lower one 
respectively) are plotted. In region (iii) {i.e., pmax > 0.31), 
Mo is taken from 0.1802 to 0.1801 in steps of 0.0001 (solid 
lines). The upper, middle and lower thin dashed curves corre- 
spond to the relation between M and pmax for d = 2.0, 1.3125 
and 1.5, respectively, (b) n = 1.25 case. From left to right, 
the Mo of each curve changes from 0.2164 to 0.2174 in steps 
of 0.0002 (solid lines) in region (i) {i.e., pmax < 0.14). In re- 
gion (ii), two long dashed lines with Mo = 0.2175 and 0.2176 
(upper one and lower one respectively) are plotted. In region 
(iii) {i.e., pmax > 0.14), Mo is taken from 0.2174 to 0.2170 in 
steps of 0.0002 (solid lines). The upper and lower thin dashed 
curves correspond to the relations between M and pmax for 
d = 1.875 and 1.3125, respectively. 
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